Comparing different coarse-grained potentials for star polymers 
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We compare different coarse-grained single-blob models for star polymers. We find that phe- 
nomenological models inspired by the Daoud-Cotton theory reproduce quite poorly the thermody- 
namics of these systems, even if the potential is assumed to be density dependent, as done in the 
analysis of experimental results. Using the numerically determined coarse-grained potential, we also 
determine the minumum value fc of the functionality of the star polymer for which a fluid-solid 
transition occurs. By applying the Hansen- Verlet criterion we find 35 < /c < 40. This result 
is confirmed by an analysis that uses the modified (reference) hypernetted chain method and is 
qualitatively consistent with previous work. 
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I. INTRODUCTION 

Soft materials are physical systems of great interest 
because of their many applications. Star polymers, ob- 
tained by tethering / polymer chains to a central micro- 
scopic core, represent a very interesting example. Indeed, 
by changing the functionality /, one can interpolate be- 
tween linear chains (/ — 1,2) and colloidal spheres, cor- 
responding to / ^ 1. Moreover, their particular nature 
is responsible for the appearance of many static and dy- 
namic peculiar features which are not observed in sus- 
pensions of hard-sphere colloids or of linear chains.-i^ 
Full-monomer simulations of star polymers are very dif- 
ficult because of the complexity of the structure of these 
molecules and of the large number of monomers involved. 
Thus, coarse-grained (CG) models, in which each poly- 
mer is represented by a single monoatomic molecule — 
sometimes it is called blob — have been widely usedi^ To 
obtain a faithful representation of the star-polymer solu- 
tion thermodynamics one should consider n-body inter- 
actions among the CG molecules^^i^ A considerable sim- 
plification occurs if one only considers the dilute regime 
in which the polymer packing fraction $ = AwRgN/ {3V) 
(here Rg is the zero-density radius of gyration, while N is 
the number of star polymers in the volume V) is at most 
of order 1. Indeed, in this case, overlaps are rare and 
the contributions due to the n-body potentials are small. 
Hence, a reasonable approximation is obtained by only 
considering the pair potential between the CG molecules, 
which is defined as 



(1) 



where the average is over all pairs of isolated star poly- 
mers such that the distance between their centers is 
r; Uint is the total intermolecular energy. Standard 
renormalization-group arguments indicate that, if one 
considers stars made of fL monomers, the adimensional 
potential /3V(r;f) is a universal function of 6 = r/Rg 



in the limit L oo. Universality implies that, for any 
given /, the limiting function /3V(6; /) is independent 
of the microscopic model. The first studies of star poly- 
mer solutions at finite density based on CG models^ used 
phenomenological potentials that were inspired by the 
Daoud-Cotton model^S They diverged logarithmically as 
6 0, as predicted theoretically^ and showed an ex- 
ponential (Yukawa) decrease for b ^ oo. Such a large- 
distance behavior is consistent with the Daoud-Cotton 
model predictions and was somewhat confirmed by the 
results of Ref. which found their experimental data 
to be consistent with an effective interaction decaying 
as e~°'^. Such a behavior was also confirmed — albeit 
with quite large errors — by numerical simulationsii of 
systems with / < 50. Subsequent numerical work — 
but again arms were quite short — did not confirm the 
Yukawa behavior for / < 18^^ This is not surprising, 
since the Daoud-Cotton model^ does not apply for small 
values of /. Hence, the phenomenological potential of 
Ref. 7 cannot be used for star polymers with a small 
number of arms. For these reasons, Ref. [l^ suggested 
that the potential of Ref. should only be used for 
/ > 10. For smaller values of /, a second phenomenolog- 
ical potential was proposed. ^3. It has the correct logarith- 
mic short-distance behavior and shows a Gaussian large- 
distance decay, in agreement with the renormalization- 
group predictions obtained for linear chains (/ = 2)ji^ A 
direct numerical determination of the universal pair po- 
tential /3V(5; /) was undertaken by Hsu and Grassberger 
(HG).^^ By means of a large-scale simulation of an opti- 
mal model (the Domb- Joyce modeli^^ at a specific value 
of the interaction parameter) they obtained accurate es- 
timates of the pair potential for several values of / in the 
range 2 < / < 35 and provided an accurate parametriza- 
tion of their results which satisfied all theoretical con- 
straints. They checked the predicted logarithmic diver- 
gence for small distances^ and found a purely Gaussian 
large-distance decay for all values of / investigated. In 
particular, while at short distances their numerical po- 
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tential was close to the potential proposed in Ref. [ij, 
significant differences were observed at large distances, 
as a consequence of the different (Yukawa vs Gaussian) 
decay. 

The phenomenological potentials introduced in Refs.0, 
I13I have been extensively used to study the phase diagram 
of star polymers, ^^^^^ binary star polymer systemSj^ 
mixtures of star polymers and colloidsSl and of star and 
linear polymersj^i^ and structural arrest in dense star 
polymer systemsi^i^ Since they are quite different from 
the exact one derived by HG^^ one may question the 
quantitative and/or qualitative validity of the results ob- 
tained. It is thus worthwhile to repeat these calculations 
by using all CG models, comparing the results obtained 
by using the phenomenological potentials with those ob- 
tained by using the accurate expression of the pair po- 
tential obtained by HG,^- which we take as reference po- 
tential. 

In this paper, we study the thermodynamic behav- 
ior of dilute star-polymer solutions, by using the CG 
model based on the HG accurate pair potential. We 
determine the first virial coefficients, the pressure in 
the dilute regime, and the intermolecular structure fac- 
tor which is needed to compare the theoretical results 
with the experimental ones obtained in scattering ex- 
periments. These results are then compared with the 
analogous ones obtained by using the phenomenological 
potentials of Refs. [tUTsI. Finally, we investigate the phase 
diagram of star polymer solutions, identifying the range 
of values of / for which a liquid-solid transition occurs, 
again using the HG pair potential. 

The paper is organized as follows. In Sec. |lT] we in- 
troduce the three different CG models. In Sec. IIIII we 
discuss the thermodynamical behavior: first, we compute 
the second and third virial coefficient for each CG model, 
then the compressibility factor and the intermolecular 
structure factor. In Sec. IIII Dl we also discuss an exten- 
sion of the model of Ref. JSj, which uses a density depen- 
dent corona diameter. Finally, in Sec. II VI we discuss the 
phase diagram and in SeclVlwe present our conclusions. 



II. THE EFFECTIVE PAIR POTENTIALS: 
DEFINITIONS 

The universal pair potential /3V(6;/) was determined 
numerically in Ref. [T^ for several values of / between 2 
and 35. The final results were parametrized as 



PVHcib-J) = - In 



Vwp{b)^bfln{af/b) 



Vcib) = CfC 



-dfb^ 



(2) 



(3) 



The potential depends on five constants a/, 6/, c/, df, 
Tf which are reported in Ref. O For 6-^0, V{b; /) 
diverges?s2£ as In 1 /b with a prefactor that can be ex- 
pressed in terms of the partition- function exponents 7/, 



which are known with good precision.-'' Parametrization 
^ satisfies this property, fiV HQ{b; f) ~ 6/ln(a//6), the 
coefficient 6/ being equal to the theoretically predicted 
value. For 6 — >■ 00 the potential behaves as Cfe~'^f'' , 
where df varies between 0.405 (/ = 2) and 0.68 (/ = 35). 

Numerical studies of the properties of star polymers 
have often relied on phenomenological expressions for the 
pair potential. In Ref. 7 the following potential was pro- 
posed: 
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R<1, 
i?>l, (4) 



where Kf = 1/(1 + ^/J/2), R = r/a, and a is the 
so-called corona diameter. For i? — > the potential 
shows the expected logarithmic behavior 6/ In 1/i? with 
bf — 5/'^/2/i8_ The coefficient bf can be compared with 
the theoretical result 6/ obtained by using the accurate 
estimates of the partition- function exponents 7/1^ We 
obtain bf = 2.42(1), 9.90(3), 57.3(6) for / = 4,10,30, to 
be compared with bf — 2.22, 8.78, 45.6 for the same val- 
ues of /. Differences increase with / and range from 8% 
for / = 4 to 20% for / — 30. They are, however, expected 
to be largely irrelevant for the thermodynamics in the di- 
lute regime, in which overlaps are rare. For i? — )• 00 the 
potential behaves as e"^^^"^ / R, which is quite different 
from the behavior observed in numerical simulations (at 
least for / < 35), see Eq. ©. 

Potential ^ is expected to be reliable only for large 
values of /, for / > 10, say, i.e. in the regime to which 
the Daoud-Cotton model* applies. For small values of / 
a different potential was postulated: ^'^ 



PV2{R;f) 
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18 V ' 2^ 

5J3/2 



36r2 



-T^(i?,^-l) 



i? < 1, 

i? > 1, (5) 



where again R = r/a. The adimensional parameter r 
determines the large-r behavior of the potential and was 
determined only for f — 2 and / = 5: r = 1.03 and 
T = 1.12 in the two cases. 

In the following we will compare the predictions of the 
three CG models defined above: model MHG based on 
the exact pair potential ^ and models Ml and M2 based 
on potentials (jlj and ([5]), respectively. 



III. RESULTS: STRUCTURE AND 
THERMODYNAMICS 

A. Zero-density results 

Knowledge of the pair potential V{b; f) allows us to 
compute the universal combination A2 = B2/R^, where 
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TABLE I: In the second, third, and fourth column we report 
virial coefficient ratios computed by using potential Vhg(^; /) 
(HG). In the last three columns we report literature values 
obtained from full-monomer (FM) simulations: a refers to 
Ref. m, h to Ref. IH, c to Ref. HO, d to Ref. [Sjie to Ref. [H. 
For a review of older estimates of A2, see Ref. l3a. 



f 


A2(HG) ^3(HG) 5(HG) 


A2(FM) 


A3(FM) g(¥M) 


2 


5.51 


4.97 


0.164 


5.500(3)" 


9.80(2)" 0.324(1)" 


4 


10.04 


31.8 


0.316 


9.979(9)'' 


39.56(16)'' 0.397(2)' 


5 


12.25 


54.5 


0.363 






6 


14.66 


85.6 


0.398 


14.174(16)'' 90.1(0.4)'' 0.449(2)'' 


10 


20.62 


204 


0.479 






12 


23.97 


288 


0.501 


23.5(2)" 




18 


29.60 


473 


0.540 


29.6(2)^^ 


0.547(5)^* 












0.58'= 


24 


36.06 


721 


0.557 


34.21(10)^* 


661(6)^* 0.564(4)'' 


30 


38.32 


831 


0.565 


37.65(6)^* 


813(7)'' 0.574(6)'' 



the second virial coefficient B2 is defined by the expan- 
sion of the (osmotic) pressure 11, 

-4^ = l + B2P + i33p' + 0(p3), (6) 

in powers of the concentration p = N/V. Indeed, A2 is 
related to the pair potential by the exact relation 

A2{f)^2i, dbb'(l-e-^''^'-^f^). (7) 

In TableUvife report the estimates of A2 for several values 
of / obtained by using potential If we compare these 
results with those obtained in the literature from differ- 
ent full-monomer simulations [column A2(FM)] we ob- 
serve reasonable agreement (the largest deviation, 6.5%, 
is observed for / = 24), confirming the adequacy of 
parametrization ([2]) for all values of /. 

Using model MHG we can also compute the universal 
combinations ^3 = B^R~^ and g = B3 / B2 = A3/A2, in- 
volving the third virial coefficient B3. Since three-body 
interactions are neglected in the CG model, these esti- 
mates differ from those that would be obtained in the 
exact, full-monomer polymer model. Therefore, the ob- 
served discrepancies give us quantitative indications of 
the role of the neglected many-body interactions in CG 
single-blob star-polymer models. The results reported in 
Table [T] show that many-body forces apparently become 
less relevant as / increases. For instance, if we consider 
the relative deviation Ag = 1 — g(IIG)/.g(FM) for the g 
parameter, we find a very large discrepancy for f — 2, 
Ag « 50%, but only Ag « 11% for / = 6. For / > 18, 
the CG model apparently reproduces the full-monomer 
results. This is particularly encouraging since it implies 
that CG models provide increasingly better approxima- 
tions as / increases, i.e. exactly in the regime in which 
full-monomer simulations become unfeasible. 



TABLE II: Estimates oia/Rg, A2, and A3 for models Ml and 
M2. 







Ml 




M2 


f 

J 


a/Rg{Ml) AaiMl) g{Ml) 


a/Rg{M2) A3{M2) g{M2) 


z 


0.807 


2.01 


0.067 


1.333 D.OO U./09 


4 


0.894 


16.4 


0.163 




5 


0.933 


30.6 


0.204 


1.349 61.9 0.412 


6 


0.974 


50.9 


0.237 




10 


1.062 


141 


0.332 




18 


1.195 


376 


0.429 




30 


1.328 


722 


0.492 
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FIG. 1: Comparison of the CG potentials of models MHG, 
Ml, and M2: (top) / = 5; (bot tom) / = 30. For / = 30 we 
fix r = 1 in model M2 (see Sec HITU)) . 



Let us now perform the same analysis for models Ml 
and M2. Since the potentials Vi(r;/) and V2{r; f) de- 
pend on the corona diameter a, to obtain quantitative 
predictions we must determine the ratio cr/Rg. As al- 
ready suggested in Ref.l2ll we fixa/Rg so that all poten- 
tials give the correct result for the second- virial combina- 
tion A2. This guarantees that models Ml and M2 have 
the correct thermodynamic behavior as $ — >■ 0. Using 
the estimates of A2 reported in Table U [column yl2(IIG)] 
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we obtain the estimates oi a/Rg reported in Table HIl For 
/ = 18, Ref. obtained a/Rg ~ 1.26 from the analysis 
of the experimental data, while numerical simulationsii 
indicate that a/Rg ~ 1.3 is a good approximation for 
model Ml at least for / large. Our numerical results are 
fully consistent with this approximation. As long as / is 
larger than 10, the expected range of validity of model 
Ml, a/Rg ~ 1.3 holds with a relative error of at most 
20%. Such an approximation also holds for model M2 
for the two values of / we consider. In the following, we 
consider models Ml and M2 using the values of a/Rg 
reported in Table |lll Hence, by construction, all mod- 
els (Ml, M2, and MHG) have the same thermodynamic 
behavior for $ — > 0. In Fig. [1] we compare the differ- 
ent potentials for / = 5 and / = 30. For / = 5 poten- 
tials M2 and MHG are very close in the whole interesting 
range h = r/Rg < 2.5. For larger values, potential M2 
decreases slightly faster: PVhg decays as 1.76e~°'^'^''^, 
while (3Vm2 decays as 4.34e~"-^^'' . Significant differences 
are instead observed for potential Ml, both for small val- 
ues of 6 — it underestimates PVhg — and for large values 
of 6, where it decays slower. Similar discrepancies are 
observed for / = 30. 

Let us now compare the third-virial coefficient combi- 
nations ^3 and see Tables HI and HIl For model Ml, 
discrepancies are quite large for / small, confirming the 
inadequacy of the parametrization for these values of /. 
For / > 10 discrepances are smaller, but still not negli- 
gible. For / = 18 .g(Ml) differs from g(HG) and ^(FM) 
by 21%. For / = 30 the discrepancy decreases to 13%. 
For / < 10 model M2 should be used. Also in this case, 
we observe significant discrepances from the results ob- 
tained by using the MHG model, but, at least for / = 2, 
model M2 appears to provide a better approximation to 
the full- monomer results. 



B. Integral-equation methods 

In order to study the finite-density behavior we 
used integral-equation methods.'^'* As usual in these ap- 
proaches we considered the pair distribution function 
g{v), the corresponding correlation function /i(r) = 
g{r) — 1, and the direct correlation function c(r) defined 
by the Ornstein-Zernike relatior^^ 



h{v) = c(r) + p (fis c{s)h(r - s 



(8) 



This equation must be supplemented by a closure rela- 
tion. Three different closures were used to check the ac- 
curacy of the results. We first used the hypernetted chain 
(HNC) equation34 



ff(r) 



-c(r) 



(9) 



which is known to be quite accurate for soft interactions. 
In our case potentials diverge as r — > 0, hence a better 
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FIG. 2: Relative deviations AZhnc ~ Zhnc /Zry — 1 and 
AZmhnc ~ Zmhnc/Zry — 1 for / = 30 and model MHG. 
The virial route (HNC,v and MHNC,v) and the compressibil- 
ity route (HNC,c and MHNC,c) have been used. 



approximation should be provided by the Rogers- Young 
(RY) closurei^s 



5(r) = e-^^("-) 1 + 



1 



fir) 



,(h(r)-c(r))/(r) 



1 



(10) 



with 



f{r) = l-e-"\ (11) 

The consistency parameter a was redetermined at each 
density by requiring the equality of the compressibility 
computed by using g{r) (compressibility route) and that 
computed by using the virial pressure. This was done 
itcratively until the relative difference between the two 
quantities was less than 0.1%. We found that aRg in- 
creases with varying between 0.5 and 2-2.5 as $ in- 
creases from 0.3 to 1. 

For / > 30 — in this case the potential has a quite 
hard core — we also used the reference or modified HNC 
(MHNC) method. ^'''^^ The closure relation is written as 

5(r) = e-'^^W+''W-^M+^W, (12) 

where E{r) is the bridge function. For E{r) we used the 
bridge function of a system of hard spheres of diameter 
d at the same density (it can be computed quite pre- 
cisely by using the results reported in Refs. 37.38). The 
diameter d, or equivalently the hard-sphere packing frac- 
tion TjHS = 7T(fip/6, was determined by using the Lado 
criterion)^ which is a way to implement thermodynamic 
consistency between the virial and the energy route: 

dE{r;riHs) 



d-^r[g{r) - gHs{r;riHs)]- 



0. 



(13) 



driHS 

To check the accuracy of the results we computed the 
compressibility factor 

P 



Z 



ksTp 



(14) 
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for model MHG by using the HNC and the RY closure; 
for / > 30 also the MHNC closure was used. As an ex- 
ample, in Fig. [2] we show t^Z^f^c = Zhnc /Zry — 1 and 
^Zmhnc = Zmhnc /Zry — 1, where Zhnc, Zry, and 
Zmhnc are the compressibility factors computed by us- 
ing the three different methods. Here we take / = 30 
and use potential MHG. For the HNC closure, we ob- 
serve a significant difference between the pressure com- 
puted by using the virial and the compressibility route, 
which however decreases as $ increases. For the MHNC 
closure, results are much more consistent, the difference 
being at most 5%, again decreasing as $ becomes large. 
The HNC and MHNC virial pressure differ only slightly 
from the RY result. If the compressibility route is used, 
differences are larger. They however decrease as the den- 
sity increases. The somewhat large difference between 
the two thermodynamic routes observed for the HNC 
method is due to the somewhat hard core of the poten- 
tial, that diverges as r — >■ 0. It is smaller for smaller 
values of /: for / = 6 and 18, Zhnc,c and Zhnc,v differ 
at most by 2.2% and 14%, respectively. In all cases the 
RY closure appears to be reliable, with an error that is 
probably of the order of a few percent and that decreases 
as $ increases. 

The results we will present in the following have been 
obtained by using the RY closure. We will use again 
the MHNC closure in Sec. IIV[ where we will discuss the 
fluid-solid transition in star-polymer solutions. 



C. Comparing the potentials at finite density 
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FIG. 3: Relative deviation AZ = Z/Zhg - 1 for model Ml 
(/ = 10, 18, 30) (top) and model M2 (/ = 2, 5, 18, 30) (bot- 
tom) as a function of the polymer volume fraction <&. Zhg is 
the compressibility factor obtained by using potential ([2]). 



Let us now discuss the behavior of the different mod- 
els at finite density in the dilute regime $ < 1. Let 
us first consider the compressibility factor Z. In Fig. [3] 
we report the relative deviations AZ = Z/Zhg ~ 1: 
where Z is computed in models Ml and M2, and Zhg 
is computed in model MHG. As should be expected on 
the basis of the results for A3 and g, model Ml under- 
stimates the true CG compressibility factor Zhg- For 
$ = 1, we find deviations AZ = -24%, -31%, -32% 
for / — 10, 18, 30, which are quite significant. Model M2 
instead overestimates Zhg- Differences are, however, sig- 
nificantly smaller: for $ = 1, we find AZ = 2%, 6% for 
/ = 2,5. 

Since model M2 has the observed Gaussian decay and 
appears to be relatively accurate, we tried to check 
whether it is possible to extend it to other values of 
/, beside / = 2, 5. For this purpose we should fix 
both ajRg and t. If we require g(M2) sa g(HG) be- 
side y4.2(M2) « A2(HG), we obtain r « 1 for all values 
of / in the range 10 < / < 35. The ratio cr/i?g is always 
consistent with 1.2. More precisely, for r = 1, we have 
ajRg = 1.203,1.202 for / = 18,30, respectively. The 
corresponding M2 potential has a large-distance behav- 
ior which is consistent with that of potential Indeed, 
we have ^V2(i?;/) ~ e'^-^'^^'' for large b = r/Rg to be 
compared with l3VHG{b; f) ~ e-'^f''\ dj = 0.65,0.68 for 



/ = 18, 30, respectively. Potential f3V2{R', f) is compared 
with pVnGib', f) in Fig. [T] On the scale of the figure, 
no significant differences are observed, confirming that 
model M2 is a better approximation to the star-polymer 
CG potential than model Ml. In Fig. [3] we also report 
AZ for model M2 with t = 1 and / = 18,30. Again 
model M2 overestimates Zhg 1 but differences are only of 
order 10% at most. 



D. Effective potentials with density-dependent 
corona diameter 

Potential ^ was validated by comparing the theo- 
retical predictions with experimental data for the total 
scattering intensity I{q)- In the dilute regime, I{q) can 
be factorized'*^ as I{q) « P{q)S{q), where P{q) is the 
single-polymer form factor and S{q) is the intermolecular 
center-of-mass structure factor, which can be computed 
by using the three different CG models. In the CG model, 
S{q) can be obtained as ^(q) = 1 -I- ph{q), where h{q) is 
the Fourier transform of the correlation function, which 
is obtained directly in the integral-equation calculation. 

Results for $ = 1 are reported in Fig. U] for / — 
2, 5, 18, 30. The results obtained by using model M2 are 
quite similar to those obtained by using model MHG, 
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FIG. 4: Structure factor as a function of qRg for $ = 1: a) / = 2 (models M2, MHG); b) / = 5 (models M2, MHG); c) / = 18 
(models Ml, M2, MHG); d) / = 30 (models Ml, M2, MHG); For / = 18, 30 the results for model M2 are obtained by setting 
T = 1. 



both for / = 2,5, which belong to the original validity 
range of the potential, and for / — 18, 30 (we take t = 1 
in this case). On the other hand, results for model Ml 
differ significantly. The positions of the minima and max- 
ima are the same for all potentials — not surprising since 
they are fixed by the dimensions of the polymer — but po- 
tential Vi gives rise to significantly stronger oscillations 
for / > 18. 

The results presented here for model Ml are in ap- 
parent contradiction with the existing literature. In- 
deed, potential (|H) has been extensively used to ana- 
lyze experimental data, finding in all cases very good 
agreement. However, it should be noted that in all 
these comparisons a density-dependent corona diameter 
is assumed. For instance, in the original paper— deal- 
ing with 18-arm polyisoprene in methylcyclohexane, a is 
fixed by (t(<I>) = 1.26i?g(<f>), where i?g($) is the density- 
dependent radius of gyration. A careful study of the 
density dependence of the corona diameter is presented 
in Ref. 41 using 57-arm polybutadiene. They find (see the 
inset of their Fig. 1) that cr($) decreases as $ increases 
and that it behaves as for $ > 1, in agreement 

with the density scaling predicted by the Daoud-Cotton 
model- Potential (|4]) has also been validated by using 
block-copolymer micelles with / = 63.^^ However, as in 
previous cases, a density-dependent corona diameter is 
assumed. 
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FIG. 5: Log-log plot of the density-dependent corona diam- 
eter a{^)/Rg for / = 18 and / = 30. The solid line corre- 
sponds to a behavior 



To understand the role of a density-dependent corona 
diameter, we repeat the same analysis as done in exper- 
imental work. We take the MHG structure factor as the 
reference one and determine a density-dependent cr($) 
such that the Ml structure factor reproduces the MHG 
one at the same value of Since the main difference 
between the MHG and Ml structure factors is the heigth 
of the peaks, see Fig. 01 we fix (t($) by requiring that 
the peak for the lowest value of q has the same height 
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FIG. 6: Compressibility factor Z as & function of $ (left) and structure factors as a function of qRg for $ = 1 (right). Data for 
model Ml with the density-dependent corona diameter reported in Fig. [5] and for model MHG. Top: / — 18; bottom: / = 30. 



in both models. The results for / = 18 and / = 30 
are reported in Fig. [5l In both cases, cr(<I>) increases for 
$ < 0.5, then shows a maximum at a value which is 
significantly larger than that obtained by matching the 
second virial coefficients, and finally decreases. The be- 
havior looks similar to that observed in Ref. |4l| (see their 
Fig. 1). It is also roughly consistent with the predictions 
of the Daoud-Cotton model, which predicts an interme- 
diate regime in which cr($) ~ $~3/4 i-,efQj-g ^j^g onset of 
the large-$ regime where cr($) ~ <I>~^/^. Model Ml with 
corona diameter (t($) well reproduces the MHG struc- 
ture factor, as can be seen from Fig. [G] 

The computation of the pressure in the presence of 
a density-dependent potential requires particular care4^ 
We use the HNC closure (thermodynamic consistency 
does not hold, hence we cannot implement the RY ap- 
proach) and the compressibility route (the virial route 
does not provide the correct result, see Ref. H^)- The 
results are shown in Fig. [51 It is evident that Z for 
model Ml is quite different from that of model MHG. 
But, even worse, the Ml predicted compressibility factor 
shows an unphysical decrease for $ > 1, a consequence 
of the quite rapid decrease of <j{^)- Clearly, use of a 
density-dependent corona diameter worsens the thermo- 
dynamic behavior of the model. 

It is not surprising that S{q) is well reproduced while 
large differences are observed for Z. Indeed, the pair dis- 
tribution function and, therefore, also S{q) are not very 



sensitive to the large-distance behavior of the potential: 
as discussed, for instance, in Ref. [H, visibly different po- 
tentials may produce structures with essentially identical 
pair distribution functions. On the other hand, thermo- 
dynamic quantities are very sensitive to the tail of the 
potential, hence differ for models MHG and Ml, even 
when a density-dependent corona diameter is used. 

IV. STAR POLYMER PHASE DIAGRAM 

Since star polymers interpolate between linear chains, 
which only have a fluid phase for all densities, and hard 
colloids, which have a fluid-solid transition, star polymers 
are expected to behave in both ways depending on f^idl 
For f < fc one expects only a fluid phase, while fo r / > fc 
a fluid and a solid phase are expected. In Refs. fc 
was estimated to be of order 100. In Ref. [H a much more 
careful analysis was performed using model Ml, finding 
fc ~ 34. For fc < f ^ they found a small range 
of densities (note that $ « 3.4?7, where rj is the volume 
fraction defined in Ref. [Tsh ■ 1.5 < <& < 2, in which a solid 
bcc phase occurs, with reentrant melting as increases. 
For larger values of /, the solid phase was more complex, 
with several crystalline states appearing at different den- 
sities. The presence of a solid phase for large values of / 
was later confirmed experimentally: Ref. 19 observed a 
solid bcc phase by using / = 67 starlike block copolymer 
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FIG. 7: Structure factor maximum Smax. = maxS{q) as a FIG. 8: Effective hard-sphere packing fraction riHS deter- 
function of $ for / = 30, 35, 40. We use the MHG model and mined by using the MHNC closure as a function of $. The 
the RY closure. horizontal line corresponds to the boundary of the fluid-solid 

coexistence rjHS ~ 0.49. 



micelles. 

Since the thermodynamical properties of model Ml 
are quite different from those of model MHG, we wish 
now to check if and how the conclusions of Ref. [l^ 
change when potential ([2]) is used. To determine the 
presence of a liquid-solid transition we use two different 
approximate methods. First, we use the Hansen- Verlet 
criterion:'*^''*^ the phase transition occurs when the max- 
imum ^max of S{q) exceeds 2.85. Second, we use the 
MHNC method'^^'^- and identify the phase transition by 
requiring rjHS ~ 0.49, where rjHS is the packing fraction 
of the reference hard-sphere model. 

Since potential MHG has only been computed up to 
/ = 35 and the region where the fluid-solid transition oc- 
curs corresponds to / > 30, we have extended the MHG 
potential to / = 40. For this purpose we need to fix the 
five parameters that appear in Eq. Parameter 6/ is 
knowufi^ bf = 94.6. Then, we note that af and df have 
a tiny dependence on /, hence it should be safe to use a 
simple extrapolation of the values appropriate to / = 30 
and 35: a/ = 1.78, df = 0.68. Parameter r/ is not pre- 
cisely known and is approximately constant for / > 10. 
We take r/ — 0.5. To fix c/, we require the model to 
reproduce the second-virial combination A2: A2 ~ 41.96 
for / = 40.31 This gives c/ = 84.2. 

In Fig. [7] we report the maximum S'max of S{q) for 
/ = 30,35,40. For $ < 5, S'max shows a nonmonotonic 
behavior with two maxima, one in the dilute region and 
one for $ ^ 1. For / = 30 and 35 the maximum S'max is 
always smaller than 2.85, hence no fluid-solid transition 
is expected. For $ = 40 instead, Smax is larger than 
2.85 in the density range 1.6 ^ ^ 2.4 (the maximum 
corresponds to <I> « 1.94 with S'max = 3.13). Hence, the 
Hansen- Verlet method allows us to infer that 35 < /c < 
40. Moreover, the solid phase should appear at values of 
<& close to 2, i.e. in a range of densities that is similar to 
that reported in Ref. [H. 

The same analysis can be repeated by using the MHNC 
closure. In Fig. [5] we report rjns as a function of $ for 



/ = 30,35,40. For / = 30,35 the effective hard-sphere 
packing fraction is always lower than 0.49, which gives 
the boundary of the fluid-solid coexistence line. For / — 
40, instead it reaches a maximum rjus — 0.51 for $ = 
1.81. Hence, the MHNC analysis predicts crystallization 
in a small $ interval that extends from ^ — 1.^ to ^ — 
2.0. 

The results of the two analyses are fully consistent and 
allow us to conclude with confidence that 35 < /c < 40. 
We therefore confirm the conclusions of Ref. [l^ although 
we predict the fluid-solid transition to occur for slightly 
larger values of / (Ref. [H predicted /c ~ 34). 



V. CONCLUSIONS 

In this paper we investigate the thermodynamic be- 
havior of three different CG models appropriate to de- 
scribe dilute star-polymer solutions. Model MHG uses 
the exact pair potential, model Ml, which should be 
applied for / > 10, is inspired by the Daoud-Cotton 
model, while model M2 is a phenomenological modifi- 
cation which shows a large-distance Gaussian behavior 
and which is expected to be realistic for / < 10. We 
find that model M2 provides a reasonable approxima- 
tion to the thermodynamics and to the structure, even 
for / > 10. On the other hand, model Ml significantly 
underestimates the pressure and does not provide the 
correct structure factor. If one wishes to reconcile struc- 
tural results for model Ml with those of model MHG, one 
might consider a density dependent corona diameter, as 
it is usually done in the analysis of the experimental data. 
However, while good agreement is obtained for the struc- 
ture factor, model Ml with state-dependent interactions 
is completely inconsistent from a thermodynamic point 
of view: the compressibily factor shows an unphysical 
decrease as the density increases. 

We investigate in detail the phase diagram of star poly- 
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mers by using model MHG. We use both the Hansen- 
Verlet criterion and the MHNC approximation to esti- 
mate the smallest value fc of the functionality / for which 
a fluid-solid transition occurs. Both analyses are consis- 
tent with 35 < /c 40. Our findings are in qualitative 
agreement with those of Ref. [3 which predicted fc ~ 34 
by using model Ml, and with experiments, which ob- 



served crystallization for / > 60. Our data are also con- 
sistent with the presence of reentrant melting for / close 
to 40: the solid phase should be stable only in a small 
density interval centered around $ « 1.8-1.9. 

We thank Giuseppe D'Adamo for useful comments. 
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